#include "TRandom.h"
void HistoCluster()
{
gROOT->Reset();
#include "Riostream.h"
#include <string>
#include <map>

const unsigned int Nevts = 20000;
const double min_point_energy = 5.;//eV per deposition point
const double min_protein_energy = 5.;//eV in one protein to damage it
const double minE_damage_RNA = 5.;//eV perdeposition point in RNA
const double maxE_damage_RNA = 5.;// eV it is used for the linear damage function with the min value
const int Nmax_proteins = 250;
const double RNA_GeomXS = 0.03;
const double Nucleocapsid_GeomXS = 0.50;

Int_t evt,copy, evt_temp=0;
Double_t x, y, z, e;
Double_t damagedprotein_perevt=0;

TCanvas *c1 = new TCanvas("Iron 25 GeV", "Iron 25 GeV",200,10,700,500);
c1->Divide(2,2);
c1->SetLogy();
//---------------------------------------------
//-------Spikes Protein damage-----------------
ifstream in;
in.open("_P.out");
TH1D *Ph = new TH1D("","", Nmax_proteins, 0, Nmax_proteins);
TH1D *PDamaged_perevt = new TH1D("","", 50, 0, 50); //Protein Max
int c= 0;
while(1){
	in>>evt>>x>>y>>z>>e>>copy;
	if(!in.good()) break;
	if(evt!=evt_temp)
		{
		int damagedproteins=0;
		for(int i=0; i<Nmax_proteins; i++)
			{
			if(Ph->GetBinContent(i)>min_protein_energy)damagedproteins++;
			}
		PDamaged_perevt->Fill(damagedproteins);
		Ph->Reset();
		evt_temp=evt;	
		}

	if(e>min_point_energy)Ph->Fill(copy,e);
	}
	PDamaged_perevt->SetYTitle("Probability");
	PDamaged_perevt->SetXTitle("Number of damaged Proteins per incident primary");
	PDamaged_perevt->SetTitle("Spike Protein damage");
	Double_t ProteinDamageMean = PDamaged_perevt->GetMean();
	c1->cd(1);
	gPad->SetLogy();
	PDamaged_perevt->Scale(1./PDamaged_perevt->Integral());
	PDamaged_perevt->Draw();
//Ph->Scale(141./Nevts);
//Ph->Draw();


//------------------------------------------------------------
//-------Membrane Glycoprotein Protein damage-----------------
ifstream inMP;
inMP.open("_MP.out");
TH1D *MPh = new TH1D("","", Nmax_proteins, 0, Nmax_proteins);
TH1D *MPDamaged_perevt = new TH1D("","", 20, 0, 20);
c= 0;
while(1){
	inMP>>evt>>x>>y>>z>>e>>copy;
	if(!inMP.good()) break;
	if(evt!=evt_temp)
		{
		int damagedproteins=0;
		for(int i=0; i<Nmax_proteins; i++)
			{
			if(MPh->GetBinContent(i)>min_protein_energy)damagedproteins++;
			}
		MPDamaged_perevt->Fill(damagedproteins);
		MPh->Reset();
		evt_temp=evt;	
		}

	if(e>min_point_energy)MPh->Fill(copy,e);
	}
	MPDamaged_perevt->SetYTitle("Probability");
	MPDamaged_perevt->SetXTitle("Number of damaged Proteins per incident primary");
	MPDamaged_perevt->SetTitle("Membrane GlycoProtein damage");
	Double_t MProteinDamageMean = MPDamaged_perevt->GetMean();
	c1->cd(2);
	gPad->SetLogy();
	MPDamaged_perevt->Scale(1./MPDamaged_perevt->Integral());
	MPDamaged_perevt->Draw();



//------------------------------------
//--------RNA-------------------------
const int maxRNAdamage_perevt = 100;//set by user upon estimation always higher than protein max
const int maxCapsiddamage_perevt = 300;
ifstream in1;
in1.open("_M4.out");
c= 0;
TRandom *R = new TRandom(time(0));
TH1D *RNAh = new TH1D("","", maxRNAdamage_perevt, 0, maxRNAdamage_perevt);
TH1D *Ncaph = new TH1D("","", maxCapsiddamage_perevt, 0, maxCapsiddamage_perevt);
evt_temp=0;
int RNAdamageN = 0;
int NcapsiddamageN =0;
while(1){
	in1>>evt>>x>>y>>z>>e;
	
	if(!in1.good())break;
	if(evt!=evt_temp)
		{
		RNAh->Fill(RNAdamageN);	
		Ncaph->Fill(NcapsiddamageN);	
		RNAdamageN=0;
		NcapsiddamageN=0;
		evt_temp=evt;
		}
	Int_t damage=0;
	if(e>minE_damage_RNA)
	{
	if(e>maxE_damage_RNA)damage=1;
	else 	{
		Double_t slope = 1./(maxE_damage_RNA-minE_damage_RNA);		
		Double_t intercept = -slope*minE_damage_RNA;
		Double_t randam = R->Rndm();
		if(randam < (slope*e+intercept))damage=1;
		}

	Double_t value = R->Rndm();
	if(value <= RNA_GeomXS){if(damage)RNAdamageN++;}
	else if(value <= (RNA_GeomXS+Nucleocapsid_GeomXS))NcapsiddamageN++;
	}
	
	}


	RNAh->SetYTitle("Probability");
	RNAh->SetXTitle("Number of RNA damage per incident primary");
	RNAh->SetTitle("RNA damage");
	Double_t RNADamageMean = RNAh->GetMean();
	c1->cd(4);
	gPad->SetLogy();
	RNAh->Scale(1./RNAh->Integral());
	RNAh->Draw();	

	Ncaph->SetYTitle("Probability");
	Ncaph->SetXTitle("Number of NCapsid damage per incident primary");
	Ncaph->SetTitle("NCapsid damage");
	Double_t NCapsidDamageMean = Ncaph->GetMean();
	c1->cd(3);
	gPad->SetLogy();
	Ncaph->Scale(1./Ncaph->Integral());
	Ncaph->Draw();	

double ssRNA=0.;
double ssProt=0.;
double ssMProt=0.;
for(int i=2; i<maxRNAdamage_perevt; i++)
{
ssRNA+=RNAh->GetBinContent(i);
ssProt+=PDamaged_perevt->GetBinContent(i);
ssMProt+=MPDamaged_perevt->GetBinContent(i);
}
cout<<"0000000000000000000000000000000000000000000000000000000000000"<<endl;
cout<<"-------------------------------------------------------------"<<endl;
cout<<"-------------------------- RNA Damage -----------------------"<<endl;
cout<<"No RNA damage "<<RNAh->GetBinContent(1)<<endl;
cout<<"One RNA damage and above 1(included)  "<<ssRNA<<endl;
cout<<"Average RNA damage  "<<RNADamageMean<<endl<<endl;

cout<<"---------------------------------------------------------------"<<endl;
cout<<"------------------------------ Averages -----------------------"<<endl;
cout<<"Average spike Proteins damaged  "<<ProteinDamageMean<<endl;
cout<<"Average number of damaged Sproteins per RNA break   "<<ProteinDamageMean/RNADamageMean<<endl<<endl;
cout<<"Average Membrane GlycoProteins damaged  "<<MProteinDamageMean<<endl;
cout<<"Average number of damaged Mproteins per RNA break   "<<MProteinDamageMean/RNADamageMean<<endl<<endl;
cout<<"---------------------------------------------------------------"<<endl;

/*
cout<<"----------------------------Integrals------------------------"<<endl;
cout<<"Total Proteins damage probability/incident primary "<<ssProt<<"  Total RNAdamageProb/incident primary  "<<ssRNA<<endl;
cout<<"Ratio total protein damage / RNA damage  "<<ssProt/ssRNA<<endl<<endl;
*/


}

